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Abstract 

In (Hansen et al. 2002) we presented a new approach for measuring non-Gaussianity of the 
Cosmic Microwave Background (CMB) anisotropy pattern, based on the multivariate empirical 
distribution function of the spherical harmonics a£ m of a CMB map. The present paper builds 
upon the same ideas and proposes several improvements and extensions. More precisely, we exploit 
the additional information on the random phases of the a£ m to provide further tests based on the 
empirical distribution function. Also we take advantage of the effect of rotations in improving the 
power of our procedures. The suggested tests are implemented on physically motivated models of 
non-Gaussian fields; Monte-Carlo simulations suggest that this approach may be very promising 
in the analysis of non-Gaussianity generated by non-standard models of inflation. We address 
also some experimentally meaningful situations, such as the presence of instrumental noise and a 
galactic cut in the map. 
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I. INTRODUCTION 



Much information is expected from ongoing and future experiments aimed at measuring 
the CMB radiation. In this framework, a considerable amount of attention has been drawn 
in the recent literature by Gaussianity tests on CMB data. Such tests have been motivated 
by the prediction of the so-called theory of inflation, currently the most popular theory of 
the very early universe. Standard inflation predicts that the temperature fluctuations of 
the CMB, should be (very close to) Gaussian distributed (see the reviews On the 

other hand more recent, non-standard models of inflation predict small non-Gaussianities; 
detecting such effects would be important for understanding the physics of the primordial 
epochs [J , 0, 0, O, ■ More radical deviations from Gaussianity are predicted by models 
where topological defects are induced by phase transitions (See Ref. 

[9] 

for a review). 

Detecting non-Gaussianities is also relevant as a tool for controlling systematic effects in 
CMB data. Among the many procedures advocated in the l iterature, many are based upon 
topological properties of spherical Gaussian fields [l 



on harmonic space properties 
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22j For empirical results, see 



In a recent paper |27( we proposed a new procedure to detect non-Gaussianity in harmonic 



space using so-called empirical processes 



28, 



More precisely, let T(9, <p) denote the 



CMB fluctuations field, which we assume, as always, to be homogeneous and isotropic, for 
< 6 < 7t, < tp < 2tt. Assuming that T(9, (p) has zero mean and finite second moments, 
it is well-known that the following spectral representation holds: 



T(9, if) = 2J a imYim(Q, 



1=1 m=-i 

where Ye m (9,<p) denotes the spherical harmonics. The random coefficients (amplitudes) 
{a£ m } have zero-mean with variance (|a^ m | 2 ) = Cg. They are uncorrelated over i and |m|: 
(a£ m a* e , m ,) = CfSfS™, and a^ m = a} _ m . The sequence {Cg} denotes the angular power 
spectrum of the random field and the asterisk complex conjugation. Furthermore, if T(9, tp) 
is Gaussian, the {aim} have a complex Gaussian distribution. Upon observing T(9, (p) on 
the full sky, the random coefficients can be obtained through the inversion formula 



/7T /*7T 
/ T (9, <p)Y; m (9,<p) sin 9d9d(p , m = 0,±1,...,±Z , / = 1,2,... . (1) 
•7T JO 



Our idea in that paper was to study the empirical distribution function for the {ai m } , 
and to use these results to implement tests for non-Gaussianity in harmonic space, in the 
presence of an unknown angular power spectrum. Although the Monte Carlo evidence 
produced in (^tJ was rather encouraging, several important issues were left for further 
research. Indeed, it must be noted that although the joint distribution of the spherical 
harmonic coefficients is invariant under rotations (under Gaussianity), their numerical value 



is not; as a result, the sample outcome of the test in [27[ is not rotationally invariant. In 
this paper this remark is exploited to implement an extended, more powerful version of the 
empirical process procedure, based upon the combination of results from several different 



choices of coordinate systems. Furt 
in four different setups: whereas in 



rermore, empirical process ideas are also implemented 



2?J we focused only on the squared spherical harmonic 
coefficients, here we consider also the real and imaginary parts as separate random 
variables, which allows us to gain further efficiency improvements. More important, we 
consider the empirical process on the random phases of the Gaussian field. This approach 
is especially convenient, because under Gaussianity not only the random phases are model 
independent (compare 0|) but the Smirnov uniformization (see below) is exact, i.e. 
no bias arises. Moreover the joint information on the random phases and the random 
amplitudes uniquely identifies the distribution of the field, under Gaussianity. We address 
also some experimentally meaningful situations, such as the presence of instrumental noise 
and a galactic cut in the map; concerning the latter, we advocate a solution which may 
have some independent interest for other harmonic space methods of statistical inference. 
Our tests are implemented on models of non-Gaussianity based upon ongoing research on 
non-standard inflationary models. 

The plan of the paper is as follows: in section |ll]) we review the empirical process 
method; the extensions provided in this paper are presented in (|III|) . Section ()IV|) presents 
the results from Monte-Carlo simulations on 100 realizations of a non-standard inflationary 
model, with varying levels of non-Gaussianity. Comments and directions for further research 
are collected in the final section (TV*!) . 
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II. THE EMPIRICAL PROCESS METHOD 



The details of the empirical process approach to detect non-Gaussianity in the CMB 



were given in |27|J. In short, the method consists of a family of tests which focus on the 
total distribution of ci£ m and checks for dependencies between k £-rows. The first step is to 
transform the spherical harmonic coefficients into variables ue m which have an approximate 
uniform distribution between and 1, given that the ai m were initially Gaussian distributed. 
This is done using the Smirnov transformation, defined as 

urn = $1 (^^-] , uem = $2 f 2 ^H ^ , m = 1, 2, I, 1 = 1, 2, ...L , 



Ci J V Ci 

where $ n is the cumulative distribution function of a x 2 with n degrees of freedom and 
Ci are the power spectrum coefficients estimated from the data. The error introduced by 
using estimated Cp instead of the real underlying Cg is dealt with using a bias-subtraction, 
as described in |2j|. 

Then the joint empirical distribution function for row I is formed, 

1 1 \ k 1 

Fi.d+A^ioci, ...,a k ) = 22 ) - ai ) nUut+At,i-i,m+A mi <<Xi)> , A mi > , 

v ' m=0 I i=2 J 

where A^, determines the spacing between the rows for which the dependencies are tested 
and A mi denotes the difference in m for row i. The parameters ol^ run over the interval [0, 1]. 
The empirical process is expressed using the centered and rescaled i^...^+A ffc _ 1 given as 

r * 

Gt.i+Ai^iai, ...,a k ) = y/(£ + l) < F^j+^^ax, ...,a k ) ~ \\ a t 

I i=i 

The intuition behind this procedure is as follows: if the a£ m s are Gaussian, G converges to 
a well-defined limiting process, whose distribution can be readily tabulated. On the other 
hand, for non-Gaussian a^ m s {F(a>i, a k ) — Yli=i a i) an( i thereby G will take 'high' values 
over some parts of a-space. Thus, the analysis of some appropriate functional of G can be 
used to detect non-Gaussianity. To combine the information over all multipoles into one 
statistic, we define 

1 [(L-Ak-Jr) ^ 

K L (a u ...,a k ,r) = —= ^2 G W+A<*_i(ai, • •■,«*) , 

V L, - lAe^k-i e=1 



where L is the highest multipole where the data is signal dominated. 

The method can then be summarized as follows: the distribution of sup|A"k| is found 
using Monte-Carlo simulations of Gaussian distributed ae m . Then, for a given observed 
set of ag m , the value k max = sup|-?Tl| is found and compared to the distribution obtained 
from Monte-Carlo. The consistency of the data with a Gaussian distribution can then be 
estimated to any suitable cr-level. 

III. IMPROVEMENTS OF THE EMPIRICAL PROCESS METHOD 

The process described in the previous section can be improved in several ways. One 
important thing to notice is that if the a^ m are non-Gaussian distributed, their joint distri- 
bution will change under a rotation of the sky, and therefore it will not be invariant with 
respect to the choice of coordinate system. In spherical harmonic space, the rotation of the 
sky with the Euler angles ($1, 0, $2) can be written as, 



where the rotation matrices L>^ TO /($i, 0, $2) are described in Appendix (jXJ). This new sets 
of ai m will yield a different (although not independent) value of /c max ; averaging the /c max 
obtained from the empirical process over several rotations may thus increase the probability 
for detecting non-Gaussianity in the data. In Table HI we indeed prove that this effect is 
noticeable, considering a non-Gaussian inflationary model with Jnl = 300. To obtain this 
data, we used 100 Gaussian Monte-Carlo simulations to find the distribution of /c max - Then 
we used 100 realizations of the non-Gaussian model and for each realization we averaged 
the fc max value obtained over 45 rotation in the ©-direction. The results are shown using 
the process to test Gaussianity and dependency between 1, 2 and 3 £-rows. We define the 
lcr, 2cr and 3a detection levels as the fc max values over which we had 32%, 5% and 1% of 
the hits in the Gaussian simulations respectively. 



As seen from equation (j2J), a rotation in the ^-direction will only change the phase of 
the ai m for which the empirical process as defined above is insensitive. The reason is that 




(2) 
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TABLE I: INFLATIONARY MODEL, f NL = 300, test A 



number of rows: 


1 


2 


3 


45 ©-rotations: 


1 


2 


3 


la 


28% 


35% 


44% 




32% 


41% 


52% 


2a 


5% 


7% 


14% 




4% 


22% 


23% 


3a 


3% 


4% 


6% 




1% 


8% 


10% 



we define the ui m using the modulus of the ae m . But by defining a complex m m for which 
the real and imaginary parts are just the uniformized real and imaginary parts of the a^ m , 
the phase information can be kept, and the test will be sensitive also to rotations in the $ 
direction. The new ue m can be obtained by a similar Smirnov transformation, 




= 1,2,...,/, Z = 1,2,...L, 



where G is the cumulative distribution function of a Gaussian variable with zero mean and 
unit variance and the r and i indices are used to indicate real and imaginary parts. 



With these new ui m the standard empirical process can be applied, now also allowing the 
check for dependencies between real and imaginary parts. We will in this paper arrange the 
U£ m in two ways. First we will double the number of rows, putting the imaginary rows in 
between the real rows, then we apply the standard empirical process using Al = Am = 250 
for L = 500 as in j^J (from now on we will label this the C test). Otherwise we will just 
double the length of the £-rows, putting the imaginary parts on the negative m side. Then 
we apply the standard empirical process still with A£ = Am = 250 for L = 500 (called the 
B test). In both ways we check dependencies between real and imaginary parts. In Table ILT1 
we show results from these two new tests applied to 100 realizations of the /nl = 300 model 
used above (to be compared with the results of test A without rotations in Table |1J). Even if 
the number of detections in each test vary, the specific realizations detected by each of the 
tests also vary. This suggests that by combining the results from all tests, the number of 
detections can be increased. We suggest to make the average /c max over all 3 tests as a new 
statistic. The average has to be weighted with the mean fc max for each test, obtained from 
Monte-Carlo. Thus the new statistic would be: 



L.A i.A 
iA |_ mm i B , mm i C 
"'max ' max ' Tq max' 



TABLE II: INFLATIONARY MODEL, f NL = 300, test B,C and A+B+C 



rows(B):: 


1 


2 


3 


rows(C): 


1 


2 


3 


rows(A+B+C): 


1 


2 


3 


la 


31% 


32% 


41% 




30% 


31% 


38% 




31% 


38% 


50% 


2a 


8% 


8% 


13% 




4% 


7% 


12% 




6% 


6% 


17% 


3a 


1% 


2% 


2% 




1% 


1% 


2% 




1% 


1% 


2% 



where the letters indicate test A, B and C. 



As we showed above, the old test can be improved significantly using rotations in the 
0-directions. The two new tests will also be sensitive to rotations in the ^-directions, as 
these rotations only change the phases of the a^ m . In Table ITTTl we show a similar increase 
in the number of detections using 30 rotations in the ^-rotation. We also show the total 
number of detections, combining test A averaged over 30 0-rotations and test B and C 
averaged over 30 ^-rotations. Clearly the highest number of detections can be obtained 
combining all 3 tests with as many rotations as possible for each, and in this way testing 
all possible dependencies in spherical harmonic space. 

Finally the same sort of ideas can be implemented by working directly on the random 
phases by taking the uniformized coefficients to be ue m = arc tan(a^, m /a^ m )/7r. From a 
theoretical point of view, this choice is particularly appealing because the uniformization 
is exact. Random phases as a test for non-Gaussianity were considered also by 0|; these 
authors however do not consider the empirical process approach. We label this test D and 
we report the related results in Table IVTT1 

In practice, a relevant issue is the coupling between spherical harmonic coefficients in- 
troduced by partial sky coverage from a galactic cut. It turns out that this can be well 
controlled if the Monte-Carlo calibration of the fc max distribution are done with the same 
galactic cut. The introduction of the galactic cut has the effect of reducing the number of 
detections, but this can be significantly improved if the gap is refilled, using the two borders 
between the gap and the observed area as a mirror. This creates a discontinuity at the 
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TABLE III: INFLATIONARY MODEL, f NL = 300, test B,C and A+B+C with 30 rotations 



rows(B):: 


1 


2 


3 


rows(C): 


1 


2 


3 


rows(A+B+C): 


1 


2 


3 


la 


33% 


31% 


67% 




28% 


30% 


50% 




30% 


42% 


75% 


2a 


8% 


5% 


29% 




6% 


4% 


18% 




8% 


15% 


32% 


3a 


2% 


3% 


15% 




1% 


3% 


16% 




4% 


10% 


17% 



galactic equator, but improves the power of the procedure. In our test, using this refilling 
technique, we found no noticeable reduction in the number of detections with and without 
a galactic cut of 5° x 2 in the 100 realizations used in the simulations. 



IV. DETECTION OF NON-GAUSSIAN INFLATIONARY MODELS 

In this section we will show the results of the 4 versions of the empirical process method 
applied to 100 realizations of a non-standard inflationary model. Indeed, as discussed in 
8J, one way of parameterizing the possible presence of non-Gaussianity in the primordial 
gravitational potential $ is to represent it in the following way, 

$ = + /Wl(0 2 -<0 2 >), (3) 

where is a zero mean Gaussian random field and f^L is a non-linear parameter which can 
be observationally constrained. Our empirical examples are based on this model, using for 
simplicity a pure Sachs- Wolfe angular power spectrum (Q oc £(£ + 1)). It should be noted 
that because <fi is of the order 10 -5 , the quadratic term is five orders of magnitudes smaller 
than the linear one. 

We will use 3 levels of non-Gaussianity, f NL = 100, /jvl = 300 and Jnl = 1000 (a map 
for each J'nl with its corresponding histogram is shown in figures ^EJ)- On these realizations 
we put non-uniform noise corresponding to the noise-level as expected from the Planck LFI 
100GHz channel. We also include a 5° x 2 galactic cut. The Monte-Carlo for determination 
of the distribution of fc max will when not otherwise indicated be made using the same maps 
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MAP 




-a.aae-Ol ^^^^^^^^^^^^^^f 3.S9e-01 mK 



FIG. 1: Map with a non-Gaussian part /jvx = 100. 



MAP 




-3.76e-01 ^^^^^^^^^^^^^H- a.76e-01 mK 



FIG. 2: Map with a non-Gaussian part /nl = 300. 

with f NL = 0. 

The results of the tests are reported in Tables II VII VI II We see that the 3-rows version 
of the test provides the best performance. Among the 4 tests, none appears to detect 
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MAP 




FIG. 3: Map with a non-Gaussian part f^L = 1000. 




FIG. 4: Histogram of the map at /jvl = 100. The shaded curve is a Gaussian fit. 

significantly more than others. Down to a level /jvl = 300, more than half of the maps are 
detected at the 2a level. However at Jnl = 100, the detection level at 2a has dropped to 
about 30% at the best. Nevertheless, as the spectrum we used went only to multipole t = 500 
and that radiative transfer was not taken into account, we expect that these numbers can 
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FIG. 5: Histogram of the map at f^L = 300. The shaded curve is a Gaussian fit. 
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FIG. 6: Histogram of the map at /jvl = 1000. The shaded curve is a Gaussian fit. 



improve significantly. 




TABLE IV: TEST B 



Inl 




100 






300 






1000 




rows 


1 


2 


3 


1 


2 


3 


1 


2 


3 


1 a 


42% 


34% 


47% 


32% 


37% 




83% 


28% 


70% 


100% 


2 a 


8% 


7% 


13% 


6% 


10% 




59% 


10% 


38% 


98% 


3 a 


2% 


3% 


9% 


1% 


4% 




52% 


1% 


30% 


98% 


TABLE V: TEST C 


Inl 




100 






300 






1000 




rows 


1 


2 


3 


1 


2 


3 


1 


2 


3 


1 a 


32% 


33% 


57% 


36% 


43% 




93% 


38% 


84% 


100% 


2 a 


11% 


8% 


12% 


9% 


15% 




68% 


22% 


78% 


99% 


3 a 


2% 


2% 


3% 


3% 


12% 




45% 


7% 


69% 


99% 



V. COMMENTS AND CONCLUSIONS 

In j27| we argued that the empirical process method enjoyed a number of advantages 
over existing methods, such as being completely model independent, allowing for a rigorous 
statistical analysis, detecting the location of non-Gaussianity in harmonic space, being easily 
amendable to extensions and generalizations. In this paper the potentials of this approach 



TABLE VI: TEST B+C 



fNL 




100 






300 






1000 




rows 


1 


2 


3 


1 


2 


3 


1 


2 


3 


1 a 


40% 


34% 


50% 


36% 


41% 


95% 


40% 


87% 


100% 


2 a 


5% 


8% 


29% 


6% 


16% 


87% 


13% 


77% 


99% 


3 a 


2% 


6% 


10% 


1% 


11% 


68% 


2% 


70% 


99% 
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TABLE VII: TEST D 



Inl 




100 






300 






1000 




rows 


1 


2 


3 


1 


2 


3 


1 


2 


3 


1 a 


35% 


37% 


45% 


39% 


40% 


83% 


44% 


47% 


97% 


2 a 


9% 


4% 


19% 


16% 


5% 


72% 


30% 


17% 


97% 


3 a 


0% 


0% 


7% 


2% 


1% 


57% 


9% 


11% 


97% 



are further investigated. Three further classes of tests in the same spirit as in the previ- 
ous paper are proposed. Changes of coordinates are implemented to make the procedure 
rotationally invariant and more powerful against non-Gaussian alternatives. We have also 
considered experimentally meaningful situations such as the presence of gaps and noise. Fi- 
nally, we have tested the power of these procedures against physically reasonable models 
such as those expected by non-standard inflation; we show that for the model considered, 
significant detections of non-Gaussianities are achieved for non linear coupling parameters 
down to about 100. Overall we believe that empirical process methods emerge as valuable 
tools for further studies on the non-Gaussianity of the CMB. 
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APPENDIX A: ROTATION MATRICES 



A spherical function T(8, d>) is rotated by the operator D(a/3 , ~f) where a/37 are the three 
Euler angles for rotations |3Jj and the inverse rotation is D(— 7 — /3 — a). For the spherical 
harmonic functions, this operator takes the form, 

£ 

Y em (0'A')= ^%M7)^(M), (Al) 

m'=-£ 

where -D^, m has the form 

KfJflM) = e m ' a d e mlm (P)e^. (A2) 
Here d e m , m (j3) is a real coefficient with the following property: 

d m >m((3) = dm m i(—f3). (A3) 
The D-functions also have the following property: 

{a 2 MD e m „ m {a 1 i3 1 j 1 ), (A4) 

m" 

where (a/37) is the result of the two consecutive rotations (ai/3i7i) and (02/^272) • 
The complex conjugate of the rotation matrices can be written as 

^ram , = (~l) m+m ^(-m)(-ra')' (A5) 
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